To install Destin2, run:
devtools::install_github("yuchaojiang/Destin2/package")
This vignette will demonstrate the utility of the Destin2 package through analysis of a single-cell ATAC-seq dataset of human peripheral blood mononuclear cells (PBMCs) provided by 10x Genomics. Links to the raw data as well as code used for pre-processing and Signac-based downstream analysis can be found on the corresponding Signac vignette. For brevity, we omit the details here; we encourage the invested reader to examine the Signac vignette.
First, we load in Destin2 as well as the other packages we will be using for this vignette.
library(Destin2)
library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Hsapiens.v75)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TFBSTools)
library(JASPAR2020)
library(cisTopic)
library(mogsa)
library(ggplot2)
library(patchwork)
library(clustree)
library(slingshot)
Next, we load in the peak/cell matrix and cell metadata and perform quality control pre-processing. The following files are used in this vignette, available through 10X Genomics:
set.seed(1234)
#Reference genome
genome=BSgenome.Hsapiens.UCSC.hg19
ensdb=EnsDb.Hsapiens.v75
# load the RNA and ATAC data
counts <- Read10X_h5("../data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "../data/atac_v1_pbmc_10k_singlecell.csv",
header = TRUE,
row.names = 1
)
fragpath <- "../data/atac_v1_pbmc_10k_fragments.tsv.gz"
chrom_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = 'hg19',
fragments = fragpath,
min.cells = 10,
min.features = 200
)
pbmc <- CreateSeuratObject(
counts = chrom_assay,
assay = "peaks",
meta.data = metadata
)
# get gene annotationsg
annotations <- GetGRangesFromEnsDb(ensdb = ensdb)
seqlevelsStyle(annotations) <- "UCSC"
Annotation(pbmc) <- annotations
pbmc <- NucleosomeSignal(pbmc)
pbmc <- TSSEnrichment(pbmc, fast = F)
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
pbmc <- subset(
x = pbmc,
subset = peak_region_fragments > 3000 &
peak_region_fragments < 20000 &
pct_reads_in_peaks > 15 &
blacklist_ratio < 0.05 &
nucleosome_signal < 4 &
TSS.enrichment > 2
)
Destin2 provides methods for conducting and integrating cross-modality analyses. In order to do so, we consider the following unimodal analyses.
The first analysis method will be to normalize and reduce the peak/cell matrix using latent semantic indexing (LSI), followed by graph-based clustering.
pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc) # This by default returns an lsi reduction
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
The second unimodal analysis method will be a gene activity matrix derived from this scATAC-seq data. Here we use the GeneActivity function provided by Signac; we note that there exist alternatives such as MAESTRO that could also be used.
gene.activities <- GeneActivity(pbmc)
# This can be substituted with MAESTRO's regulatory potential model
# add the gene activity matrix to the Seurat object as a new assay and normalize it
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
object = pbmc,
assay = 'RNA',
normalization.method = 'LogNormalize',
scale.factor = median(pbmc$nCount_RNA)
)
For downstream analysis, we reduce the dimension of his matrix using PCA.
DefaultAssay(pbmc) <- 'RNA'
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
pbmc <- ScaleData(pbmc, features = rownames(pbmc))
pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc),
reduction.name = 'activity.pca', reduction.key = 'activitypca_', verbose = FALSE)
Downstream analysis oftentimes makes use of transferred cell type labels. We follow the procedure outlined in the Signac vignette.
pbmc_rna <- readRDS("../data/pbmc_10k_v3.rds")
transfer.anchors <- FindTransferAnchors(
reference = pbmc_rna,
query = pbmc,
reduction = 'cca'
)
predicted.labels <- TransferData(
anchorset = transfer.anchors,
refdata = pbmc_rna$celltype,
weight.reduction = pbmc[['lsi']],
dims = 2:30
)
pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)
pbmc <- subset(pbmc, idents = 14, invert = TRUE)
pbmc <- RenameIdents(
object = pbmc,
'0' = 'CD14 Mono',
'1' = 'CD4 Memory',
'2' = 'CD8 Effector',
'3' = 'CD4 Naive',
'4' = 'CD14 Mono',
'5' = 'DN T',
'6' = 'CD8 Naive',
'7' = 'NK CD56Dim',
'8' = 'pre-B',
'9' = 'CD16 Mono',
'10' = 'pro-B',
'11' = 'DC',
'12' = 'NK CD56bright',
'13' = 'pDC'
)
pbmc$seurat_clusters=droplevels(pbmc$seurat_clusters)
# Remove cells with low label-transfer scores
pbmc=subset(pbmc, cells=which(pbmc$prediction.score.max>0.8))
# Remove cells that are predicted extremely rarely
table(pbmc$predicted.id)
pbmc=subset(pbmc, cells=which(pbmc$predicted.id!='Dendritic cell'))
Next, we compute a cell by motif activity matrix using chromVAR. Full details of the method can be found in the corresponding paper.
# Get a list of motif position frequency matrices from the JASPAR database
pfm <- getMatrixSet(
x = JASPAR2020,
opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)
# add motif information
DefaultAssay(pbmc)='peaks'
pbmc <- AddMotifs(
object = pbmc,
genome = genome,
pfm = pfm
)
# Computing motif activities
pbmc <- RunChromVAR(
object = pbmc,
genome = genome,
new.assay.name = 'MOTIF'
)
Similarly, we reduce this cell by motif matrix using PCA.
DefaultAssay(pbmc) <- 'MOTIF'
#Using PCA dimensional reduction on motif modality and put it into Seurat Object
#Choice of 50 PCs here is arbitrary
pca_motif<-prcomp(t(pbmc@assays[["MOTIF"]]@data))$x[,1:50]
pbmc[["motif.pca"]] <- CreateDimReducObject(embeddings = pca_motif, key = "motifpca_", assay = "MOTIF")
The final unimodal analysis we will consider is topic modeling using cisTopic. Full details for the method can be found on their github. We note that the number of topics chosen as well as other computing parameters were chosen based on the computing power at hand; one may need to adjust these as needed based on individual computing capabilities.
DefaultAssay(pbmc) <- 'peaks'
temp = pbmc@assays$peaks@counts
rownames(temp)=GRangesToString(pbmc@assays$peaks@ranges, sep=c(':','-'))
cisTopicpbmc = createcisTopicObject(count.matrix = temp)
rm(temp)
cisTopicObject <- runWarpLDAModels(cisTopicpbmc, topic=c(5, 10, 20, 25, 30, 40), seed=987, nCores=6, iterations = 500, addModels=FALSE)
cisTopicObject <- selectModel(cisTopicObject, type = 'derivative')
cisTopicObject <- runUmap(cisTopicObject, target='cell', seed=123, method='Probability')
dim(cisTopicObject@selected.model$topics)
dim(cisTopicObject@selected.model$document_expects)
pbmc[['lda']] <- CreateDimReducObject(embeddings = t(cisTopicObject@selected.model$document_expects), key = "lda_", assay = "peaks")
Destin2 offers three options for cross-modality integration:
consensus PCA, multiple CCA, and weighted nearest neighbor. These are
implemented in the functions GetConsensusPCA,
GetMultiCCA, and FindMultiModalNeighbors from
the Seurat V4 package.
GetConsensusPCA and GetMultiCCA have a
similar set of main parameters:
obj: Seurat objectreduction.list: List of dimensional reduction objects
within objdims.list: Dimensions of reduction objects to use for
the joint dimension reductionreduction.name: Name of resulting joint dimension
object, to be stored in objreduction.key: Prefix for resulting joint dimension
objectassay: Name of assay used to calculate dimensional
reduction## list of dimensional reduction objects in pbmc
reductions <- list('lsi','lda','activity.pca','motif.pca')
## dimensions to use for each dimensionality reduction object; we choose 30 for this demonstration
## we drop first PC from lsi based on correlation with sequencing depth
reduction_dims <- list(2:31, 1:30, 1:30, 1:30)
pbmc <- GetConsensusPCA(pbmc, reduction.list=reductions,
dims.list = reduction_dims,
reduction.name = 'consensus.pca',
reduction.key = "consensuspca_",
assay = "peaks")
pbmc <- GetMultiCCA(pbmc, reduction.list=reductions,
dims.list = reduction_dims,
reduction.name = 'multicca.pca',
reduction.key = "multiccapca_",
assay = "peaks")
pbmc <- FindMultiModalNeighbors(pbmc, reduction.list=reductions,
dims.list = reduction_dims, verbose=FALSE)
pbmc
## An object of class Seurat
## 108317 features across 4869 samples within 3 assays
## Active assay: peaks (87561 features, 87561 variable features)
## 2 other assays present: RNA, MOTIF
## 9 dimensional reductions calculated: lsi, activity.pca, motif.pca, lda, consensus.pca, multicca.pca, consensus.umap, multicca.umap, wnn.umap
We demonstrate the utility of the cross-modality framework that Destin2 offers through graph-based clustering and non-linear dimension reduction and visualization using UMAP.
########################
#### Consensus PCA
########################
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction='consensus.pca',
reduction.name='consensus.umap', reduction.key = 'consensusumap_', verbose = FALSE)
DefaultAssay(pbmc)='peaks'
pbmc <- FindNeighbors(object = pbmc, reduction = 'consensus.pca', dims = 1:30, verbose = FALSE, graph.name=c('consensuspca_nn','consensuspca_snn'))
pbmc <- FindClusters(object = pbmc, algorithm = 1, verbose = FALSE, graph.name='consensuspca_snn')
pbmc$consensus_clusters=pbmc$seurat_clusters
p1=DimPlot(object = pbmc, group.by='consensus_clusters', label = TRUE, reduction = 'consensus.umap') + NoLegend()+ggtitle('Consensus PCA')
p2=DimPlot(object = pbmc, group.by='predicted.id', label = TRUE, reduction = 'consensus.umap') + NoLegend()+ggtitle('Consensus PCA')
p1+p2
########################
#### MultiCCA
########################
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction='multicca.pca',
reduction.name='multicca.umap', reduction.key = 'multiccaumap_', verbose = FALSE)
pbmc <- FindNeighbors(object = pbmc, reduction = 'multicca.pca', dims = 1:30, verbose = FALSE, graph.name=c('multicca_nn','multicca_snn'))
pbmc <- FindClusters(object = pbmc, algorithm = 1, verbose = FALSE, graph.name='multicca_snn')
pbmc$cc_clusters = pbmc$seurat_clusters
p1=DimPlot(object = pbmc, label = TRUE, group.by='cc_clusters',reduction = 'multicca.umap') + NoLegend()+ggtitle('MultiCCA')
p2=DimPlot(object = pbmc, group.by='predicted.id',label = TRUE, reduction = 'multicca.umap') + NoLegend() +ggtitle('MultiCCA')
p1+p2
########################
#### Weighted NN
########################
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
## 23:03:57 UMAP embedding parameters a = 0.9922 b = 1.112
## 23:03:58 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 23:04:00 Initializing from normalized Laplacian + noise (using irlba)
## 23:04:00 Commencing optimization for 500 epochs, with 150900 positive edges
## 23:04:09 Optimization finished
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3, resolution = 2, verbose = FALSE)
pbmc$wnn_clusters = pbmc$seurat_clusters
p1=DimPlot(object = pbmc, label = TRUE, group.by='wnn_clusters',reduction = 'wnn.umap') + NoLegend()+ggtitle('WNN')
p2=DimPlot(object = pbmc, group.by='predicted.id',label = TRUE, reduction = 'wnn.umap') + NoLegend() +ggtitle('WNN')
p1+p2
For further clustering inference, the clustree allows us to build clustering trees and measure the stability of clusters as resolution increases.
#########################################################
## clustree to determine the number of clusters / resolution
#########################################################
resolution.range <- seq(from = 0.2, to = 1.2, by = 0.2)
pbmc <- FindClusters(object = pbmc, algorithm = 1, verbose = FALSE, resolution = resolution.range, graph.name='consensuspca_snn')
pbmc.cpca.clustree=clustree(pbmc, prefix = 'consensuspca_snn_res.')
p1=pbmc.cpca.clustree
# The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions.
# The stability index is automatically calculated when a clustering tree is built.
pbmc.cpca.clustree=clustree(pbmc, prefix = 'consensuspca_snn_res.', node_colour = "sc3_stability")
p2=pbmc.cpca.clustree
p1+p2
resolution.range <- seq(from = 0.2, to = 1.2, by = 0.2)
pbmc <- FindClusters(object = pbmc, algorithm = 1, verbose = FALSE, resolution = resolution.range, graph.name='multicca_snn')
pbmc.mcca.clustree=clustree(pbmc, prefix = 'multicca_snn_res.')
p1=pbmc.mcca.clustree
pbmc.mcca.clustree=clustree(pbmc, prefix = 'multicca_snn_res.', node_colour = "sc3_stability")
p2=pbmc.mcca.clustree
p1+p2
resolution.range <- seq(from = 0.2, to = 1.2, by = 0.2)
pbmc <- FindClusters(object = pbmc, algorithm = 1, resolution = resolution.range, verbose = FALSE, graph.name = "wsnn")
pbmc.wnn.clustree=clustree(pbmc, prefix = 'wsnn_res.')
p1=pbmc.wnn.clustree
pbmc.wnn.clustree=clustree(pbmc, prefix = 'wsnn_res.', node_colour = "sc3_stability")
p2=pbmc.wnn.clustree
p1+p2
The final part of our analysis uses the Slingshot
package to perform trajectory inference. The GetSlingshot
function provides a convenient wrapper to this package. A
vignette on Slingshot can be found here.
GetSlingshot takes the following parameters:
obj: Seurat objectreduction.to.construct: Dimension reduction object used
to infer lineagesreduction.to.plot: Dimension reduction object used to
cluster centroids and plot lineagescluster: Vector of cluster identity for each
samplepredicted.id: String specifying (transferred) cell
labelsgenerate.plot: If true, returns Slingshot plotsCPCALineages <- GetSlingshot(pbmc, reduction.to.construct = 'consensus.pca',
reduction.to.plot='consensus.umap',
cluster='consensus_clusters',
predicted.id='predicted.id',
generate.plot = TRUE)
CPCALineages
## class: SlingshotDataSet
##
## Samples Dimensions
## 4869 2
##
## lineages: 6
## Lineage1: 2 0 1 7 3 4 8
## Lineage2: 2 0 1 7 3 5
## Lineage3: 2 0 1 7 3 6
## Lineage4: 2 0 1 7 10
## Lineage5: 2 0 1 9
## Lineage6: 2 0 1 11
##
## curves: 0
MCCALineages <- GetSlingshot(pbmc, reduction.to.construct = 'multicca.pca',
reduction.to.plot='multicca.umap',
cluster='cc_clusters',
predicted.id='predicted.id',
generate.plot = TRUE)
MCCALineages
## class: SlingshotDataSet
##
## Samples Dimensions
## 4869 2
##
## lineages: 6
## Lineage1: 8 4 2 7 0 3 9
## Lineage2: 8 4 2 7 0 3 11
## Lineage3: 8 4 2 7 0 1
## Lineage4: 8 4 2 7 10
## Lineage5: 8 4 2 5
## Lineage6: 8 4 2 6
##
## curves: 0
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] slingshot_2.4.0 TrajectoryUtils_1.4.0
## [3] SingleCellExperiment_1.18.1 SummarizedExperiment_1.26.1
## [5] MatrixGenerics_1.8.1 matrixStats_0.62.0
## [7] princurve_2.1.6 clustree_0.5.0
## [9] ggraph_2.1.0 patchwork_1.1.2
## [11] ggplot2_3.3.6 mogsa_1.30.0
## [13] cisTopic_0.3.0 JASPAR2020_0.99.10
## [15] TFBSTools_1.34.0 BSgenome.Hsapiens.UCSC.hg19_1.4.3
## [17] BSgenome_1.64.0 rtracklayer_1.56.1
## [19] Biostrings_2.64.1 XVector_0.36.0
## [21] EnsDb.Hsapiens.v75_2.99.0 ensembldb_2.20.2
## [23] AnnotationFilter_1.20.0 GenomicFeatures_1.48.4
## [25] AnnotationDbi_1.58.0 Biobase_2.56.0
## [27] GenomicRanges_1.48.0 GenomeInfoDb_1.32.4
## [29] IRanges_2.30.1 S4Vectors_0.34.0
## [31] BiocGenerics_0.42.0 sp_1.5-0
## [33] SeuratObject_4.1.2 Seurat_4.2.0
## [35] Signac_1.8.0 Destin2_0.0.0.9000
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 scattermore_0.8
## [3] R.methodsS3_1.8.2 tidyr_1.2.1
## [5] bit64_4.0.5 knitr_1.40
## [7] irlba_2.3.5.1 DelayedArray_0.22.0
## [9] R.utils_2.12.0 data.table_1.14.4
## [11] rpart_4.1.16 KEGGREST_1.36.3
## [13] RCurl_1.98-1.9 generics_0.1.3
## [15] snow_0.4-4 RhpcBLASctl_0.21-247.1
## [17] cowplot_1.1.1 RSQLite_2.2.18
## [19] RANN_2.6.1 future_1.28.0
## [21] bit_4.0.4 tzdb_0.3.0
## [23] spatstat.data_2.2-0 xml2_1.3.3
## [25] httpuv_1.6.6 assertthat_0.2.1
## [27] DirichletMultinomial_1.38.0 viridis_0.6.2
## [29] xfun_0.34 hms_1.1.2
## [31] jquerylib_0.1.4 evaluate_0.17
## [33] promises_1.2.0.1 fansi_1.0.3
## [35] restfulr_0.0.15 progress_1.2.2
## [37] caTools_1.18.2 dbplyr_2.2.1
## [39] igraph_1.3.5 DBI_1.1.3
## [41] htmlwidgets_1.5.4 spatstat.geom_2.4-0
## [43] feather_0.3.5 rsparse_0.5.1
## [45] purrr_0.3.5 ellipsis_0.3.2
## [47] backports_1.4.1 dplyr_1.0.10
## [49] annotate_1.74.0 sparseMatrixStats_1.8.0
## [51] biomaRt_2.52.0 deldir_1.0-6
## [53] vctrs_0.4.2 ROCR_1.0-11
## [55] abind_1.4-5 withr_2.5.0
## [57] cachem_1.0.6 ggforce_0.4.1
## [59] progressr_0.11.0 checkmate_2.1.0
## [61] sctransform_0.3.5 GenomicAlignments_1.32.1
## [63] prettyunits_1.1.1 goftest_1.2-3
## [65] cluster_2.1.4 lazyeval_0.2.2
## [67] seqLogo_1.62.0 crayon_1.5.2
## [69] arrow_9.0.0.2 genefilter_1.78.0
## [71] labeling_0.4.2 pkgconfig_2.0.3
## [73] tweenr_2.0.2 nlme_3.1-160
## [75] ProtGenerics_1.28.0 rlang_1.0.6
## [77] globals_0.16.1 lifecycle_1.0.3
## [79] miniUI_0.1.1.1 filelock_1.0.2
## [81] BiocFileCache_2.4.0 doSNOW_1.0.20
## [83] polyclip_1.10-4 lmtest_0.9-40
## [85] graph_1.74.0 Matrix_1.5-1
## [87] zoo_1.8-11 ggridges_0.5.4
## [89] RcisTarget_1.16.0 png_0.1-7
## [91] viridisLite_0.4.1 rjson_0.2.21
## [93] float_0.3-0 bitops_1.0-7
## [95] R.oo_1.25.0 lda_1.4.2
## [97] KernSmooth_2.23-20 DelayedMatrixStats_1.18.2
## [99] blob_1.2.3 stringr_1.4.1
## [101] parallelly_1.32.1 spatstat.random_2.2-0
## [103] readr_2.1.3 CNEr_1.32.0
## [105] scales_1.2.1 memoise_2.0.1
## [107] graphite_1.42.0 GSEABase_1.58.0
## [109] magrittr_2.0.3 plyr_1.8.7
## [111] ica_1.0-3 gplots_3.1.3
## [113] zlibbioc_1.42.0 compiler_4.2.1
## [115] BiocIO_1.6.0 RColorBrewer_1.1-3
## [117] fitdistrplus_1.1-8 Rsamtools_2.12.0
## [119] cli_3.4.1 listenv_0.8.0
## [121] pbapply_1.5-0 text2vec_0.6.2
## [123] MASS_7.3-58.1 mgcv_1.8-40
## [125] tidyselect_1.2.0 stringi_1.7.8
## [127] highr_0.9 yaml_2.3.6
## [129] ggrepel_0.9.1 grid_4.2.1
## [131] sass_0.4.2 fastmatch_1.1-3
## [133] tools_4.2.1 future.apply_1.9.1
## [135] parallel_4.2.1 rstudioapi_0.14
## [137] TFMPvalue_0.0.8 foreach_1.5.2
## [139] AUCell_1.18.1 gridExtra_2.3
## [141] farver_2.1.1 Rtsne_0.16
## [143] digest_0.6.30 rgeos_0.5-9
## [145] shiny_1.7.2 pracma_2.4.2
## [147] Rcpp_1.0.9 later_1.3.0
## [149] RcppAnnoy_0.0.19 httr_1.4.4
## [151] colorspace_2.0-3 XML_3.99-0.11
## [153] tensor_1.5 reticulate_1.26
## [155] lgr_0.4.4 splines_4.2.1
## [157] uwot_0.1.14 RcppRoll_0.3.0
## [159] spatstat.utils_3.0-1 graphlayouts_0.8.3
## [161] plotly_4.10.0 xtable_1.8-4
## [163] jsonlite_1.8.2 poweRlaw_0.70.6
## [165] tidygraph_1.2.2 corpcor_1.6.10
## [167] R6_2.5.1 mlapi_0.1.1
## [169] pillar_1.8.1 htmltools_0.5.3
## [171] mime_0.12 glue_1.6.2
## [173] fastmap_1.1.0 BiocParallel_1.30.4
## [175] codetools_0.2-18 utf8_1.2.2
## [177] lattice_0.20-45 bslib_0.4.0
## [179] spatstat.sparse_2.1-1 tibble_3.1.8
## [181] curl_4.3.3 svd_0.5.2
## [183] leiden_0.4.3 gtools_3.9.3
## [185] GO.db_3.15.0 survival_3.4-0
## [187] rmarkdown_2.17 munsell_0.5.0
## [189] GenomeInfoDbData_1.2.8 iterators_1.0.14
## [191] reshape2_1.4.4 gtable_0.3.1
## [193] spatstat.core_2.4-4